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We present an effcient Monte Carlo method for lattice glass models which are characterized 
by hard constraint conditions. The basic idea of the method is similar to that of the JV-fold 
way method. Since it is a waste of computational time trying an insertion of a particle which 
is forbidden by constraint conditions, we just try an insertion which does not conflict with the 
constraint conditions by using a list of the insertable sites. The list is made at the beginning of 
the simulation, and it is updated whenever the particle configuration is changed. We applied the 
method to a lattice glass model proposed by Biroli and Mezard. We first evaluated the efficiency 
of the method through measurements of the autocorrelation function of particle configurations. 
As a result, we found that, when the chemical potential \x is large, the relaxation time of the 
present method is about 10 3 times shorter than that of the standard Monte Carlo method. This 
result shows that the efficiency of simulations is much improved, even given the computational 
cost to keep the list of the insertable sites. We next examined how the efficiency of the replica 
exchange method concerning chemical potential is influenced by a choise of a local update 
method. The results show that the efficiency of the replica exchange method is considerably 
improved by the use of the proposed update method. For example, when JV a it e = 1024, the 
crgodic time te, which is the average round-trip time of a replica in chemical-potential space, 
with the present local update method is more than 10 2 times shorter than that with the 
standard local update method. Lastly, we calculated the density of states of the model by 
combining the present method with the Wang-Landau method. 

KEYWORDS: Monte Carlo, lattice glass model, hard constrained problem, regular random graph, close- 
packing density 



1. Introduction 

Understanding the glass transition is one of the present 
major challenges in condensed matter physics. ^ There 
has been a long-standing controversy over 50 years on 
whether the glass transition observed in experiments is 
just a dynamical transition characterized by a drastic 
slowing down of dynamics or a true static transition with 
a thermodynamic singularity. To understand the nature 
of the glass transition, numerous studies have been done 
(see reviews 2-6 and references therein). 

In the theoretical study of the glass transition, lattice 
glass models 7-11 ) have played an important role. They 
are roughly divided into two distinct classes. The first 
class is kinetically constrained models in which glassy be- 
havior is induced by dynamical constraint rules for move- 
ments of particles. The Kob- Anderson model 7 -' is a rep- 
resentative of this class. This model exhibits a dynamical 
transition with a divergence of the relaxation time even 
though it has simple equilibrium properties. The second 
class is geometrically constrained models in which pos- 
sible particle configurations are geometrically restricted 
by hard constraint conditions. This class of models was 
first introduced by Biroli and Mezard in 2002. 8 ) We here- 
after refer the model as the BM model. The BM model 
on a regular random graph has been studied by the cav- 
ity method 12 - 1 with one-step replica-symmetry-breaking 
ansatz, 8, 13, 14 ) which is conjectured to provide an exact 
solution in equilibrium for a class of mean-field like mod- 
els on sparse random graphs. As a result, it is found that 



the model exhibits a static glass transition accompanied 
by a thermodynamic anomaly. 13,14 ) 

To unveil the properties of the actual glass, studies of 
lattice glass models in finite dimensions are effective. In 
such studies, numerical simulations play an important 
role because most of analytic methods are not applicable 
to systems in finite dimensions. However, it is obvious 
that numerical simulations in glassy systems seriously 
suffer from the inherent characteristic of the glass, i.e., 
slow dynamics. Because of this slow relaxation, it is very 
difficult to investigate equilibrium properties in glassy 
systems. This problem is resolved to some extend by us- 
ing an extended ensemble method 15 ) such as the multi- 
canonical method, 16, 17 - 1 the Wang-Landau method, 18, 19 ' 
and the replica exchange method. 20 ' However, this prob- 
lem of slow equilibration has not been completely settled 
by these methods. Furthermore, we can not use these 
methods to study dynamical characteristics of the model 
because they are crucially changed by the use of these 
methods. 

To relieve the problem of slow dynamics in lattice 
glass models, we invented an efficient Monte Carlo (MC) 
method. Although this method is applicable to both the 
kinetically and geometrically constrained lattice glass 
models, here we explain the method for simulating the 
BM model in this manuscript. The basic idea is rather 
similar to that of the A^-fold way method, 21 - ) which pro- 
vides us a rejection-free MC mainly used in classical spin 
systems. Since it is a waste of computational time trying 
an insertion of a particle which is forbidden by the con- 
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straint conditions, wc just try an insertion which does 
not conflict with the constraint conditions by using a list 
of the insertable sites. The list is made at the beginning 
of the simulation, and it is updated whenever the particle 
configuration is changed. We hereafter call the method 
as the list referring MC (LRMC) method. The LRMC 
method consists of two local updates: insertion-deletion 
update and particle-hole exchange update. Owing to the 
use of the insertion list, we can drastically increase the 
acceptance ratio of these two local updates in compari- 
son with that of the standard MC simulation. Concerning 
the particle-hole exchange update, the acceptance ratio 
is unity. As a result, the relaxation time of the LRMC 
method becomes much shorter than that of the standard 
MC method, particularly at high densities. For example, 
when the chemical potential \i in the grand-canonical en- 
semble is 6.5, the relaxation time of the LRMC method 
is about 10 3 times shorter than that of the standard MC 
method. Although wc have to pay some computational 
cost to make and update the insertion list, the efficiency 
of simulations is much improved by using the method 
because of this large reduction of the relaxation time. 
The LRMC method is applicable no matter whether the 
model is defined on a sparse random graph or a usual 
lattice in finite dimensions. Another important merit of 
the LRMC method is that, as will be shown in §4.2, the 
use of this method does not change the dynamics of the 
model essentially. This means that we can use the LRMC 
method to investigate dynamical properties of the model. 

By using the LRMC method, wc performed the fol- 
lowing calculations: We first measured the fj, dependence 
of the average occupation density by the LRMC method, 
standard MC method, and replica exchange method. The 
results are shown in §4.1. It is confirmed that correct re- 
sults are obtained by the LRMC method. We also show 
in this subsection that how the probability Plm of the 
system being at a local minimum depends on the density 
and size. We next evaluated the efficiency of the LRMC 
method in comparison with the standard MC method 
through measurements of the autocorrelation function. 
The results are shown in §4.2. We then examined how the 
efficiency of the replica exchange method 20 - 1 concerning 
chemical potential is influenced by the choice of the local 
update method. To be specific, we adopted the LRMC 
and standard MC methods as a local update method, 
and compared the two results. The results are shown in 
§4.3. We see that the efficiency of the replica exchange 
method is greatly improved by the use of the LRMC 
method. For example, when iV S it e = 1024, the ergodic 
time te, which is the average round-trip time of a replica 
in chemical-potential space, with the LRMC method is 
more than 10 2 times shorter than that with the stan- 
dard MC method. This means that the efficient local 
update method is important to make extended ensemble 
methods like the replica exchange method more effec- 
tive. Lastly, we calculated the density of states of the 
BM model by combining the LRMC method with the 
Wang-Landau method. 

The outline of the paper is as follows: In §2, we define 
the BM model. In §3, we present the LRMC. The results 
on the acceptance ratio are also shown in this section. In 



§4, we show our simulation results. Section 5 is devoted 
to conclusions. Technical details for updating the list are 
described in Appendixes. 

2. Model 

In this section, we define the BM model 8 -* that is stud- 
ied in the present work. The BM model is a kind of lat- 
tice glass models. A binary variable <Ji is defined on each 
site. The variable <Ji denotes whether a site i is occu- 
pied by a particle (cr, = 1) or not (c; — 0). In this 
work, we consider the BM model defined on a regular 
random graph. Each site is connected with k neighbour- 
ing sites which arc chosen randomly from all of the sites. 
The reason why we choose the BM model on a regu- 
lar random graph as a benchmark of the LRMC method 
is that its properties are well investigated by the cavity 
method. 8 ^ 13 ' 14 ) We now emphasize that, as mentioned 
above, the LRMC method is applicable regardless of the 
geometry of the graph. A particle configuration {<Ji\ is 
restricted by hard constraints that neighbouring parti- 
cles of each particle should be less than or equal to I. 
The BM model is characterized by the two integers k 
and I. They satisfy the inequality k > I. 

To investigate the thermodynamic properties of the 
BM model, we consider to sample a state {<7i} from dis- 
tribution P{(Ji} by MC simulations. The distribution is 
given as 

P{ctJ = Z- x C{ai}W{(n}. (1) 

In this equation, Z is the partition function defined by 
Z = Tr { ai }C '{<7i}W {(Ji\ and C{<Ji} is an indicator func- 
tion which is one if {<Ti} satisfies all of the constraint 
conditions or zero otherwise. W^{t7i} is a weight of the 
particle configuration {<Jj}. For example, in the case of 
the grand-canonical ensemble, W^lcr^} is given as 

W{ai} = exp \pN{<Ti}} , (2) 

where (i is a chemical potential and N{<Ji} is the number 
of particles defined by the equation 

^W = E- , "i. ( 3 ) 

where N s n e is the number of sites. 

In the present study, we will focus on the BM model 
on a regular random graph with k = 3 and 1 = 1. All 
of numerical simulations are performed in this model. 
In the grand-canonical ensemble, the model exhibits a 
static glass transition with a one-step replica symmetry 
breaking at fi s 6.8. 13 ' 14 ' 26) 

3. LRMC Method 

3.1 Standard MC method and its drawbacks 

In this subsection, we explain a standard MC method 
and its drawbacks. The following is the procedure of 
the standard MC method with the Metropolis transition 
probability: 22 ' 

(a) Prepare an initial state. The initial state can be 
chosen arbitrarily if it satisfies the constraint con- 
ditions. 

(b) Choose a site k at random. 
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(c) Create a new state {oft from the present state {oft 
by changing the value of o^ from to 1 or vice versa. 

(d) Accept the change into {oft with the probability 

A({a,}^{a',}) = C{v',}mm(l" 



WW,}]' (4) 

If it is not accepted, unchange the state from {oft. 

(e) Return to (b) and repeat the steps (b)-(d). 

We next explain the drawback of this standard MC 
method. We hereafter consider the grand-canonical en- 
semble whose equilibrium weight is given by eq. (2). 
Then, the acceptance ratio in eq. (4) is rewritten as 

C{oft K = 0), 



A({<J,} -»• {<xft) 



exp(-M) (o-fc = 1), 



(5) 



where we have used eq. (2) and assumed that fi is not 
negative for simplicity. We also have used the fact that 
C{oft is 1 when o^ = 1 because a deletion of a particle 
never conflicts with the constraint conditions. When [i is 
large, there are few empty sites at which we can insert 
a particle. Therefore, when ctj. = 0, the trial to insert a 
particle into the site k fails in most cases. On the other 
hand, when Ofc = 1, the acceptance ratio A({oft — > {oft) 
is quite small because [i is large. As a result of the small 
acceptance ratio in step (d) , relaxation to the equilibrium 
state becomes very slow. 

3.2 Local update I: insertion- deletion update 

In order to overcome the difficulty in the standard MC 
method, we introduce two efficient local updates for the 
LRMC method in the following two subsections. 

In this subsection, we explain the first update which 
consists of insertions and deletions of a particle. We here- 
after call it insertion- deletion update. The basic idea is 
as follows: Since it is a waste of computational time try- 
ing an insertion of a particle which is forbidden by the 
constraint conditions, we just try an insertion which does 
not conflict with the constraint conditions. In order to do 
that, we make a list of the sites into which we can insert 
a particle and update the list whenever the particle con- 
figuration is changed. We choose an insertion site at ran- 
dom from the list. The acceptance ratio of the insertion 
and that of the deletion are chosen so that the detailed 
balance condition is satisfied. We determine whether we 
try an insertion or deletion of a particle with the equal 
probability. 

As mentioned before, this method is rather similar to 
that of the A-fold way method 21 -* in the sense that we 
utilize a list to improve the simulation efficiency. It is 
worth pointing out that we have to pay some computa- 
tional cost to make and update the list. It will be dis- 
cussed in detail in §4.2 whether the LRMC method is still 
effective or not even if this additional computational cost 
is taken into account. 

We now start concrete description of the insertion- 
deletion update. We assume that we have a list of the 
sites from which we can delete a particle and that of 
the sites into which we can insert a particle. The former 
list is the same as that of the occupied sites because a 



deletion of a particle never conflicts with the constraint 
conditions. The method to detect the sites on which the 
insertion list has to be updated and the method to up- 
date the insertion list are explained in appendices A 
and B, respectively. The following is the flow chart of 
the insertion-deletion update: 

(1) Choose whether we try an insertion or deletion with 
the equal probability. 

(2a) If the insertion is chosen in step (1), select an in- 
sertion site at random from the insertion list and 
accept the insertion with an acceptance ratio 



Ai({o-i} -> {oft) = min 1 



(6) 



where Aftoft is the number of the sites into which 
we can insert a particle, Aftoft is the number of 
particles defined by eq. (3), and {oft is the parti- 
cle configuration created from {oft by inserting a 
particle into the insertion site. 

(2b) If the deletion is chosen, select a deletion site at ran- 
dom from the deletion list and accept the deletion 
with an acceptance ratio 

N{ai}W{aiY 



A R ({o-i} — » {oft) = min 1 



(7) 



AftoftWftoft 

where {oft is the particle configuration created from 
{o~i} by deleting a particle from the deletion site. 

(3) Update the insertion and deletion lists when the in- 
sertion or deletion is accepted. 

(4) Return to (1) and repeat the steps (l)-(3). 

It is straightforward to show that the procedure de- 
scribed above satisfies the detailed balance condition. We 
consider two particle configurations {o~i\ and {cr-}. The 
latter configuration {oft is created from {oft by inserting 
a particle into an insertion site. The transition probabil- 
ity T({oft — > {oft) from {oft to {oft is given as 



i 



x A z ({oft {oft), (8) 



2 Aftoft 

In the right hand side of eq. (8), the first factor is the 
probability that the insertion is chosen in step (1), the 
second factor is the probability that the proper inser- 
tion site is chosen among the K{ai] insertion sites in 
step (2a), and the third factor is the probability that 
the transition into {oft is accepted in step (3). Simi- 
larly, the transition probability of the reversal process 

T ({ a i} ~* {°ft) is S ivcn as 



x A R ({<,1} -> {„,}). (9) 



The second factor in the right hand side of eq. (9) comes 
from the fact that the proper deletion site is chosen from 
the Aftoft occupied sites. The detailed balance condition 

P{oftT({oft -> {oft) = P{oftT({oft -> {aft), (10) 

can be shown from eqs. (1), (6), (7), (8), (9), and the fact 
that C{oft = C{oft = f. 
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Fig. 1. (Color online) The chemical potential fi dependence of the 
acceptance ratio Paccept for the BM model on regular random 
graphs with k = 3 and I — I. The data measured by the LRMC 
method are denoted by full symbols and those measured by the 
standard MC method are denoted by open squares. The number 
of sites N s - ltc in the standard MC method is 512. The four bold 
arrows in the figure indicate the value of n*(N s ^ tc ) defined by 
eq. (11) at which the average number of the insertable sites is 
one. The solid line is proportional to exp(— fi). 



3.3 The size and chemical-potential dependence of the 
acceptance ratio 

To roughly examine how the efficiency of simulation is 
improved by the insertion-deletion update, we measure 
the acceptance ratio P acC ept in the grand-canonical en- 
semble by both the LRMC and standard MC methods. 
To measure Paccept, we adopted a simulated annealing 
method: The chemical potential /i is gradually increased 
from to 10 in steps of Aix = 0.01. The system is kept at 
each chemical potential for 10 5 MC steps, where one MC 
step of the LRMC method is defined by iV s it trials of the 
insertion-deletion update and subsequent A S i t0 particle- 
hole exchange updates. The particle-hole exchange up- 
date will be introduced in §3.4. The first 5 x 10 4 MC steps 
are for relaxation and the subsequent 5 x 10 4 MC steps 
are for measurement. The average over random graphs is 
taken over 100 samples. During the simulation, we also 
measure the average of K{ai\ at each fi. 

In Fig. 1, the acceptance ratio P a ccept is plotted as a 
function of li. The data measured by the LRMC method 
are denoted by full symbols and those measured by the 
standard MC method are denoted by open squares. In 
the standard MC method, we adopted the Metropolis 
transition probability. We find from the figure that the 
acceptance ratio of the LRMC method is two or three 
orders of magnitude higher than that of the standard 
MC method when \i is large. We also see that the accep- 
tance ratio of the LRMC method depends on not only 
the chemical potential // but also the size A s i te . For a 
fixed fj,, it increases with the size. In contrast, we have 
checked that the acceptance ratio of the standard MC 
method hardly depends on the size (not shown in the 
figure). 

The four arrows in the figure indicate the value of 
chemical potential li* (iV S j tc ) at which the average number 
of the insertable sites is unity. To be specific, ^*(A srte ) 
is defined by 



Fig. 2. (Color online) The acceptance ratio Paccept of the LRMC 
method shown in Fig. 1 is plotted as a function of fi — fi* . The 
inset shows the size dependence of fi* which is defined by eq. (11). 



In this equation, the overline 777 denotes the average over 
different realizations of regular random graphs and the 
bracket (■ • • ) denotes the average in the grand-canonical 
ensemble. The two subscripts of the bracket denote the 
size and chemical potential. We see from Fig. 1 that 
Paccept decays exponentially above Li*(N s it e ). We also no- 
tice that Paccept is mostly determined by the difference 
/it — fx* (iVgite) ■ In Fig. 2, Paccept is plotted as a function of 
jU — /i*(A r s itc). All of the data nicely collapse into a single 
curve. The inset of Fig. 2 shows the size dependence of 
/x*(jV s ito)- We see that fx*(N s i te ) is approximately given 
as 



M*(A si tc) « log(iV 8itc ) + C. 



(12) 



To understand this curious behavior of the acceptance 
ratio, it is useful to consider how the number of particles 
changes in the standard MC method. The number of par- 
ticles increases when an insertable site is chosen and the 
trial to insert a particle into the site is accepted. On the 
other hand, the number of particles decreases when an 
occupied site is chosen and the trial to delete a particle 
from the site is accepted. Therefore, we obtain 

AN{a i} = 

iVsite 

N{cii} 
N site 

where AN{ai} is the average of the change in the number 
of particles, A(N —> N + 1) is the probability that a trial 
to insert a particle is accepted, and A(N —> N — 1) is the 
probability that a trial to delete a particle is accepted. If 
we consider the grand-canonical ensemble and adopt the 
Metropolis acceptance ratio, eq. (13) is rewritten as 



(jV{o-<} -> N{o-i} + l) 
A(N{ai}^N{<Ti}-lj, (13) 



AN{o~i} = A - „ e 



As 



A si 



(14) 



where we have assumed that /i is positive. Since the aver- 
age number of particles does not change once the system 
is equilibrated, we obtain 

<AJV{ £ r i » M = J-({K{<Ti})» 



As 



o. 



1. 



(11) 



(15) 

On the other hand, in the grand-canonical ensemble, the 
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acceptance ratios of the LRMC method, i.e., eqs. (6) and 
(7), are given as 

Ai(W^K'})=nhn(l,^l 7r ), (16) 

AniWi} K'}) = mm (l, ^0f- ) , (17) 

where we have used eq. (2). Therefore, we notice from 
cq. (15) that the two acceptance ratios become unity 
if the number of particles N and that of the insertable 
sites K are equal to their mean values. In practice, they 
are not equal to their mean values because the system 
size is finite. However, their relative standard deviations 
become smaller and smaller as the system size increases. 
This is the reason why the acceptance ratio in Fig. 1 
increases with the size. 

We next consider the size dependence of p*. From 
eqs. (11) and (15), we obtain 

(K{ai}) Neiteili , = p(N site , fi *)N ate e-' t ' = 1, (18) 

where p = (N{ai}/N s i te ) is the average occupation den- 
sity. Since p hardly depends on the size and its p de- 
pendence is much weaker than e _M , it is appropriate 
to approximate p by a constant. Then, we easily obtain 
cq. (12). 

Lastly, we consider how the exponential decay of 
P accept above p*(N sit c) is understood. We first focus 
on the acceptance ratio for deletion. By approximating 
N{(Ji\ in cq. (17) by the average value p(N sitc , p)N sito , 
we obtain 

JV-fotfe - " w p(iV sitc ,^)iV sito e- M 

« p(N sitCl p*)N sitc e-^e-^^ 

= e~^-^\ (19) 

where we have assumed that the average occupation den- 
sity at p is close to that at p* to go from the first line 
to the second. This approximation is valid if p* is large 
enough because the average occupation density is almost 
constant for large p (see Fig. 4) and we now focus on the 
region p > p* . We have also used eq. (18) to go from the 
second line to the third. Concerning the number of the 
insertable sites K, its average value is smaller than one 
since we consider the case p > p* . It becomes smaller 
and smaller as p increases. However, because {o\-} is a 
particle configuration after a particle is removed from a 
site, if {c-} in eq. (17) is larger than one. Note that it is 
always possible to insert a particle into the deletion site. 
Therefore, when p 3> p* , K{a[} in eq. (17) is well ap- 
proximated by one. By using this fact and substituting 
eqs. (19) into eq. (17), we obtain 

A K ({cTi} -+ {oft) « e-^*\ (20) 

for p p*. We next turn to the acceptance ratio for 
the insertion. Now the point is that the average num- 
ber of particles does not change once the system is equi- 
librated. Therefore, because the insertion process and 
deletion process are chosen with the equal probability, 
the equilibrium value of and that of A\ should be the 
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Fig. 3. (Color online) Transition processes between {of"} and 
{ct?} in the particle- hole exchange update, {cj} is an intermedi- 
ate state (see text for its definition). The transition probabilities 
among the states are shown above and below the arrows. 



same. This means that A\ is also given by eq. (20). This 
is the reason why P acC ept in Fig. 1 decays exponentially 
when p > p* . On the other hand, as we discussed above, 
the acceptance ratio Paccept is close to unity when p < p* 
because K{<Ji} « (K{ai}) S> 1. From these results, one 
can naturally expect that Paccept satisfies a scaling law 

Paccept (N sitc ,p) = G[p - p*(N sito )l (21) 

where 67 is a scaling function which behaves as 

<w-{.ir j A Y J!J: <*> 

As mentioned above, the validity of the scaling is nicely 
demonstrated in Fig. 2. 

3.4 Local update II: particle-hole exchange update 

In this subsection, we explain the second local update 
which consists of exchanges of a particle for a hole. We 
hereafter call it particle-hole exchange update. The accep- 
tance ratio of this local update is unity when the weight 
W{o"j} depends only on the number of particles. The 
grand-canonical ensemble given as eq. (2) is an example 
which satisfies this condition. We hereafter consider the 
case that this condition is satisfied. As in the previous 
subsection, it is assumed that we have the insertion and 
deletion lists. The flow chart of the particle-hole exchange 
update is as follows: 

(1) Select a deletion site at random from the delctablc 
sites. Then, delete a particle from the site. 

(2) Update the insertion and deletion lists. 

(3) Select an insertion site at random from the in- 
sertable sites except the deletion site in step (1). 
Then, insert a particle into the site. 

(4) Update the insertion and deletion lists. 

(5) Return to (1) and repeat the steps (l)-(4). 

In the step (3), it is forbidden that the deletion site in 
step (1) is chosen as the insertion site so that this ex- 
change process causes a change in the particle configu- 
ration. However, in practice, this useless exchange is al- 
lowed if the deletion site in step (1) is the only insertable 
site. 

We now show that this procedure satisfies the detailed 
balance condition. We consider two particle configura- 
tions {c^} and {erf} that are transferred from each other 
by a particle-hole exchange update, as shown in Fig. 3. 
The two configurations are the same except at the two 
sites p and q. We assume that the number of particles 
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Fig. 4. (Color online) The chemical potential fi dependence of the 
average occupation density for the BM model on regular random 
graphs with k = 3 and 1 = 1. The data obtained by the LRMC 
method, standard MC method, and replica exchange method 
are denoted by open circles, open triangles, and open squares, 
respectively. The number of sites JV s ite is 512. The average over 
random graphs is taken over 100 samples. 



Fig. 5. (Color online) The probability Plm of the system being 
at a local minimum is plotted as a function of p for several sizes. 
The model is the BM model on regular random graphs with 
k = 3 and 2 = 1. The average over random graphs is taken over 
100 samples. 



and I = 1. 



in the two configurations is N'. In Fig. 3, {a\} is an in- 
termediate state which is realized after the step (1) of 
the exchange process. We see that the two sites p and 
q are empty in the intermediate state. It is important 
to notice that the intermediate state in the transition 
process from {erf} to {erf} is the same as that in the 
reverse process. As shown in Fig. 3, the transition from 
{erf} to {<rf} occurs if and only if the site q is chosen 
from the N' deletable sites in step (1) and the site p is 
chosen from the K{a\] — 1 insertable sites in step (3) 
(recall that the deletion site in step (1) is not chosen as 
the insertion site). Therefore, the transition probability 

T (K A } -> Wf}) is { N '( K Wl} - The transition 

probability of the reverse process T({af} — > {erf}) is 
also {N 1 (K{a\] — 1)} _1 in the same way. As a result, we 
find 

T{{at) -> Wf}) = T({af} -> {of}). (23) 

On the other hand, when the weight of each state de- 
pends only on the number of particles, we obtain 

PK A } = P{af}, (24) 

since C{af* } = C{af} = 1 and the number of parti- 
cles is the same in the two configurations. It is clear 
from eqs. (23) and (24) that the detail balance condi- 
tion eq. (10) is satisfied. 

It should be noted that the sampling from the distri- 
bution cq. (1) can not be achieved by only repeating the 
exchange process since it preserves the number of par- 
ticles. It is necessary to combine the insertion-deletion 
update with the particle-hole exchange update to yield 
the sampling from the distribution eq. (1). 

4. Results 

This section is devoted to show our simulation results. 
For comparison, simulation is performed not only by the 
LRMC method but also by the standard MC method. 
In standard MC simulations, we adopted the Metropolis 
transition probability. All of simulations arc performed 
for the BM model on regular random graphs with k = 3 



4-1 Results obtained by simulated annealing simulations 
To confirm that correct results are obtained by the 
LRMC method, we first measured the p dependence of 
the average occupation density by a simulated annealing 
method. The conditions of measurements are the same 
as those of acceptance ratio measurements in §3.3. Mea- 
surements are performed by both the LRMC and stan- 
dard MC methods. The number of sites A^i te is 512. The 
average over random graphs is taken over 100 samples. 

The results of measurements of the average occupa- 
tion density p(p) are shown in Fig. 4. When p is small, 
both the data coincide with each other. This shows that 
correct results are obtained by the LRMC method. On 
the other hand, for larger fj,, p(p) obtained by the LRMC 
method is clearly larger than that obtained by the stan- 
dard MC method. In Fig. 4, we also show the data ob- 
tained by the replica exchange method. The data of the 
LRMC method are slightly smaller than those of the 
replica exchange method for large p. However, their dif- 
ference is rather small. These facts indicate that the equi- 
libration is accelerated by the use of the LRMC method. 

During the simulation, we also recorded the number 
of times that a state with a density p is sampled. At the 
same time, we recorded the number of times that a local 
minimum state with a density p is sampled. The local 
minimum is here defined as a state that there is no site 
into which we can insert a particle. 23 ^ This is nothing 
but a state with K{ai} — 0. From these data, we can 
calculate the probability Plm of the system being at a 
local minimum as a function of p. We emphasize that, 
unlike the standard MC method, such measurement is 
easily feasible in the LRMC method because we always 
store the information of K{o~i\ during the simulation. 
By performing this measurement at several sizes, we ex- 
amined how this probability depends on the density and 
size. The result is shown in Fig. 5. We see that the p 
dependence of Plm becomes steeper and steeper as the 
size increases. When p is above 0.5745, Plm is almost 
unity for all of the sizes. This result implies that the 
close-packing density is around this density. 
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Fig. 6. (Color online) (i) The time dependence of the autocorre- 
lation functions C(t)'s measured by the three different methods 
and (ii) their scaling plot. The model is the BM model on regular 
random graphs with k = 3 and 1 = 1. The number of sites N s i te is 
512 and the chemical potential fi is 6.5. The average over random 
graphs is taken over 100 samples. For each sample, the average 
over thermal noise is taken over 320 independent MC runs with 
different random number sequences. In Fig. (ii), C(t)'s measured 
by the three methods (a), (b), and (c) are plotted as a function 
of t/i? a , t/Rh> ar >d t, respectively. i? a and i?b are evaluated to 
be 901 and 20.8, respectively, by fitting. 

4-2 Estimation of the efficiency 

In order to estimate the efficiency of the LRMC 
method quantitatively, we measure the autocorrelation 
function of particle configurations defined by 

c(t) = E» M* + *>i(*0)-£i M*0F m) 
E^O^-E^C*')) 2 

In the present work, the average over random graphs 
777 was taken over 100 samples and the average in the 
grand-canonical ensemble (• • • ) was taken over 320 MC 
runs with different random number sequences. Therefore, 
we performed 32000 MC runs to calculate C(t) for each 
fi. The time for equilibration, i.e., t' in eq. (25), is chosen 
to be sufficiently larger than the relaxation time of the 
autocorrelation function. 

To make comparisons between the LRMC and stan- 
dard MC methods, we measured the autocorrelation 
function by the following three methods: 

(a) Standard MC method. 

(b) Only the insertion-deletion update. 

(c) Combination of the insertion-deletion update with 
the particle-hole exchange update. 

In the method (b), one MC step is defined by iV s it e trials 
of the insertion-deletion update. In the method (c), one 
MC step is defined by iV s it e trials of the insertion-deletion 
update and subsequent iV s it e particle-hole exchange up- 
dates. The method (c) was used in the measurements of 
the acceptance ratio (Figs. 1 and 2), the average occupa- 
tion density (Fig. 4), and the probability P LM (Fig. 5). 

In Fig. 6(i), we show C(t)'s measured by the three 
methods as a function of MC steps. The number of sites 
7V S itc is 512 and the chemical potential [i is 6.5. We see 
that C(£)'s of the methods (b) and (c) decay much faster 
than that of the method (a). We next try a scaling of 
these data. The result is shown in Fig. 6 (ii) . In the fig- 
ure, C(t) of the method (c) is plotted as a function of t. 
whereas C(t) of the method (a) and that of the method 
(b) arc plotted as a function of t/R a and t/Rb, respec- 
tively. The scaling factors i? a and i?b arc evaluated by 
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Fig. 7. (Color online) The fi dependence of the scaling factor i?b 
for three sizes. The model is the BM model on regular random 
graphs with k = 3 and 1 = 1. 
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Fig. 8. (Color online) The fi dependence of the scaling factor R a 
for three sizes. The model is the BM model on regular random 
graphs with k = 3 and 1 = 1. 

fitting. We see that all of the data collapse into a single 
curve. We confirmed that this scaling holds well for all 
of the sizes (from 128 to 512) and chemical potentials 
(from to 6.5) we examined. These results indicate that 
the local updates introduced in the LRMC method do 
not change intrinsic dynamics of the system. By doing 
such analyses, we evaluated i? a and i?b for several iV s i te 
and /z. 

We see from Fig. 6(i) that C(t) of the method (c) de- 
cays faster than that of the method (b). However, since 
the computational time of the method (c) per one MC 
step is larger than that of the method (b), it is not clear 
solely from this result whether the particle-hole exchange 
update is really effective or not. Therefore, we first mea- 
sured the computational times of the two methods per 
one MC step. As a result, we found that the computa- 
tional time of the method (c) is about 2.5 times larger 
than that of the method (b) regardless of the size and 
chemical potential. This means that the particle-hole ex- 
change update is effective if the scaling factor i?b is larger 
than 2.5. Figure 7 shows how i?b depends on the size and 
chemical potential. We notice that i?b is larger than 2.5 
for all of the sizes and chemical potentials we examined. 
This shows that the particle-hole exchange update is al- 
ways effective regardless of N sit c and fi. We also find that 
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Fig. 9. (Color online) The [i dependence of the relaxation time r 
of the method (c) for three sizes. The model is the BM model on 
regular random graphs with k = 3 and I = 1. r is calculated by 
eq. (26). 
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Fig. 10. (Color online) The sample average of a maximum den- 
sity p max observed during simulation is plotted as a function of 
l/N B i te . The data obtained by the replica exchange method and 
those by the simulated annealing method are denoted by full cir- 
cles and full triangles, respectively. The straight line is a fitting 
line obtained from the data of the replica exchange method. The 
close-packing density calculated by the cavity method is indi- 
cated by the arrow. 



1) i?b increases with /z. 

2) i?b increases with decreasing N s n e . 

The fact 1) shows that the particle-hole exchange update 
is more effective when [i is large. The fact 2) indicates 
that the relative importance of the particle- hole exchange 
update increases with decreasing the size because the ac- 
ceptance ratio of the insertion-deletion update decreases 
with decreasing the size (see Fig. 1) and that of the 
particle-hole exchange update is always unity. In other 
words, the particle-hole exchange update compensates 
for the drawback of the insertion-deletion update that 
its acceptance ratio rapidly decays when \i > n*(N s it e )- 
Now let us turn to the comparison between the LRMC 
and standard MC methods. Since the method (c) is al- 
ways more efficient than the method (b), it is enough to 
compare the two methods (a) and (c). To compare the 
two methods, we first evaluated the scaling factor R a for 
several sizes and chemical potentials. Figure 8 shows the 
results. We see that the ratio increases exponentially with 
/x. This means that the superiority of the LRMC method 
to the standard MC method increases rapidly with fj,. To 



understand this behavior of i? a , it is worth recalling that 
the acceptance ratio of the standard MC method decays 
exponentially (see Fig. 1). On the other hand, such rapid 
decrease in the acceptance ratio does not exist in the 
LRMC method. As shown in Fig. 1, the acceptance ratio 
of the insertion-deletion update decays more gradually 
than that of the standard MC method. Furthermore, the 
acceptance ratio of the particle-hole exchange update is 
unity. This difference in the acceptance ratio is probably 
the main reason why the LRMC method is much more 
efficient than the standard MC method when fi is large. 

We next compared the two methods from a view point 
of computational time. As a result, we found that, if the 
standard MC method is implemented in a usual way, 
the computational time of the LRMC method per one 
MC step is about 7.0 times larger than that of the stan- 
dard MC method. The efficiency of the LRMC method 
is estimated by the number i? a divided by this ratio of 
the computational time. Therefore, Fig. 8 shows that the 
LRMC method is much more efficient than the standard 
MC method for large \i even given the smallness of the 
computational time of the standard MC method. How- 
ever, we can reduce the computational time of the stan- 
dard MC method greatly by using a multi-spin coding 
technique, 24,25 ) which is a special numerical method for 
models with discrete, especially binary, variables. Unfor- 
tunately, this technique is not applicable to the LRMC 
method. We found that, when this technique is used, the 
computational time of the LRMC method is about 230 
times larger than that of the standard MC method. This 
means that the superiority of the LRMC method might 
be reduced considerably by the use of the multi-spin 
coding technique. However, when [i = 6.5, the LRMC 
method is still about 5 times more efficient than the 
standard MC method. We also remark that the multi- 
spin coding technique is not always applicable when it is 
used with an extended ensemble method. For instance, 
the multi-spin coding technique is not incompatible with 
the Wang-Landau method that is known to be efficient 
for evaluating the density of states, because the weight 
W{oi} in eq. (2) always changes during the simula- 
tion. Meanwhile, as demonstrated in §4.4, the LRMC 
method can be efficiently coupled with the Wang-Landau 
method. 

Lastly, we examined how the relaxation time depends 
on the size and chemical potential. The result is shown 
in Fig. 9. Particle configuration {di} is updated with the 
method (c) and the relaxation time r is defined by the 
integral 



dt'C(t'), 



(26) 



where t max is the maximum time until which we mea- 
sured C(t). We confirmed that C(i max ) is almost zero 
(see Fig. 6(i)). For small n, r gradually increases with 
fi and it docs not depend on the size. In contrast, r for 
large fj, depends not only on \i but also on N s i t0 , and it 
rapidly increases with them. This result is consistent with 
a previous result obtained by the cavity method that the 
model exhibits a dynamical transition with the breaking 
of ergodicity at /id « 6.4. 13 ) This result indicates that, 
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Fig. 11. (Color online) (i) The chemical potential fj, dependence of the acceptance ratio of the replica exchange P cx for several sizes. 
The model is the BM model on regular random graphs with k = 3 and 1 = 1. The data with the standard MC method and those 
with the LRMC method are denoted by full symbols and open symbols, respectively, (ii) The chemical potential of a specific replica 
is plotted as a function of MC step. The data with the standard MC method and those with the LRMC method are shown in the 
upper and lower panels, respectively. The number of sites N a ^ te is 512. (iii) The size dependence of the ergodic time Tg. The data 
with the standard MC method and those with the LRMC method are denoted by full circles and full triangles, respectively. In the 
measurements of P ex and te, the average over random graphs is taken over 100 samples. 



in the thermodynamic limit, the relaxation time diverges 
at this chemical potential. 

4-3 Effect of local update method on the replica ex- 
change method 

As we mentioned before, extended ensemble meth- 
ods such as the multicanonical method, 16 ' 17 - ) the Wang- 
Landau method, 18,19 ^ and the exchange MC method 20 ) 
are known to be quite effective to relieve the problem 
of slow equilibration in glassy systems. The effectiveness 
of the extended ensemble methods is demonstrated in 
Fig. 10. In the figure, the sample average of a maximum 
density p ma x observed during simulation is plotted as a 
function of l/N s i te . The measurement is performed by 
both the replica exchange and simulated annealing meth- 
ods. When iVgito is small, there is no difference between 
the two data. However, we clearly see the difference for 
large A S i t0 . Pmax measured by the simulated annealing 
method is saturated to a value around 0.5742, whereas 
Pmax measured by the replica exchange method continues 
to increase with the size. Furthermore, the extrapolated 
value of p m ax in the thermodynamic limit agrees well 
with the close-packing density p s = 0.57574 obtained by 
the cavity method. 13,14 ) 

Since extended ensemble methods work well even in 
glassy systems, one may consider that the choice of lo- 
cal update method is not so important as long as we use 
an extended ensemble method. However, there is a pos- 
sibility that the efficiency of extended ensemble methods 
depends on the efficiency of local update methods. In 
fact, it has been pointed out that details of local update 
affect the efficiency of the replica exchange method. 27 ) 
We therefore examined how the efficiency of the replica 
exchange method concerning chemical potential depends 
on the choice of the local update method. The local up- 
date methods we examined are the standard MC method 
and the LRMC method (the method (c) in the previous 
subsection). In the replica exchange method, every two 
replicas at adjacent chemical potentials Hi and Hi+i are 
attempted to be exchanged per one MC step. We set the 
lowest chemical potential and the highest one at and 



8, respectively. The intermediate chemical potentials be- 
tween them are determined so that the acceptance ra- 
tio of the replica exchange P cx is roughly constant. The 
number of replicas is 64. A common set of chemical po- 
tentials is used for two simulations with different local 
update methods. To evaluate the efficiency of the replica 
exchange method, we measured the ergodic time te- This 
is defined by the average MC steps for each replica to 
move from the highest chemical potential to the lowest 
one and return to the highest one. 

In Fig. ll(i), the acceptance ratio of the replica ex- 
change P ox with the standard MC method and that with 
the LRMC method are plotted as a function of /i. We 
see that the two acceptance ratios are completely the 
same at all of the chemical potentials and sizes. We next 
show in Fig. 11 (ii) how the chemical potential of a spe- 
cific replica changes with time by the replica exchange 
process. The number of sites A s i te is 512. We clearly see 
that the movement of the replica with the standard MC 
method is different from that with the LRMC method. In 
the former case, there is a bottle-neck in the replica move- 
ment around [i « 4. In contrast, such bottle-neck does 
not exist in the latter case. Figure 11 (iii) shows the size 
dependence of the ergodic time te of the two local update 
methods. As expected from Fig. 1 1 (ii) , the ergodic time 
with the standard MC method is much larger than that 
with the LRMC method. The ratio of the former to the 
latter increases with the size and it reaches more than 10 2 
when A s ito = 1024. These results show that the efficient 
local update method is important to make extended en- 
semble methods like the replica exchange method more 
effective. 

4-4 Measurement of the density of states 

In this subsection, we show the results on the density 
of states (DOS) O(A) of the BM model. O(A) is defined 

by 

O(A') = Tr {ai} 5(N' - N{ai})C{ai}, (27) 

where C{cr,} is the indicator function which appears in 
eq. (1). To measure the DOS, we used the Wang-Landau 
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method. 18,19 ' This method is one of the standard meth- 
ods to measure the DOS. The essence of the Wang- 
Landau method is to adjust the weight Wwh(N{ai}), 
which only depends on the number of particles, so that 
P(N) becomes constant, where P(N) is the probability 
that a particle configuration with a particle number N 
is sampled. Then, because 

P(N) cx W WL (N)Q(N), (28) 

the resultant Wwl(N) is related to H(N) as 

n(N) =CW wh (N)- 1 , (29) 

where G is a proportionality constant. We will explain in 
the next paragraph how to determine it. 

The DOS is measured in the range of < p = 
N/N sitc < 0.572. The calculation of the DOS is paral- 
lelized by dividing the range into three ranges and cal- 
culating the DOS in each range independently. To be 
specific, we set the three ranges as follows: 

(a) < p < 0.55. 

(b) 0.54 < p < 0.565. 

(c) 0.56 < p < 0.572. 

It should be noted that there is overlap between two ad- 
jacent ranges. We next explain how we determine propor- 
tionality constants of the DOS of each range. We here- 
after denote them C fJl (p = a,b,c). Firstly, C a is deter- 
mined from the condition that 0(0) = 1. This condition 
comes from the fact that the number of particles is zero 
if and only if all of the sites are empty. Secondly, Cf, is 
determined so as to minimize the sum 

AS ab =Y /k (Sa(k)-S b (k)) 2 , (30) 

where S^k) is the entropy defined by 

S^k)=\og{^(k)} (p = a,b), (31) 

fi^fc) is the DOS of the range (p) and the sum runs over 
the overlap between the two ranges. C c is determined in 
the same way. 

To adjust Wwl(-^V) so that P{N) becomes flat, we 
performed a standard procedure of the Wang-Landau 
method. 18,19 ' We first introduce a function G\vl(A0 
which is related to the weight Wwl(N) as 

W wh (N) = cx P [-Gwl(A0]. (32) 

The function Gyjh(N) is modified as 

G W l{N) = Gwl(N) + AG, (33) 

after each trial to insert or delete a particle in the 
insertion-deletion update. This modification is done no 
matter whether the trial is accepted or not. We set the 
initial AG to 0.1. To check whether P(N) is flat or 
not, we accumulate the histogram H(N), which is added 
by one whenever the number of particles becomes N. 
The histogram is checked every 1000 MC steps. We re- 
gard the histogram as flat when H(N) for all N is not 
less than 80% of the mean value of H(N). If the his- 
togram is flat, AG is halved and H(N) is reinitialized 
to zero. We stop our simulation after we halve AG 20 
times. Since we set the initial AG to 0.1, the final AG is 
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Fig. 12. (Color online) The density p dependence of the entropy 
per site s(p) = log{C(p)} / JV a ite f° r five sizes. The model is the 
BM model on regular random graphs with k = 3 and 1 = 1. The 
average over random graphs is taken over 100 samples. 

0.1 x 2~ 20 » 9.5 x 10~ 8 . We checked that G wh (N) con- 
verges well in later stages of simulations. We find from 
eqs. (29) and (32) that the resultant Gwl(^V) is related 
with the entropy S(N) = \og{fl(N)} as 

S(N)=G Wh (N) + C, (34) 

where C is a constant. 

In Fig. 12, the entropy per site s(p) = S(p)/N s i tc is 
plotted as a function of p for several N s a c . We see that 
all of the data almost completely collapse into a single 
curve. s(p) starts to drop around p ps 0.32 and it becomes 
close to zero at the highest density p = 0.572 of our 
calculation. From these data, we can calculate various 
thermodynamic quantities such as internal energy, heat 
capacity, free energy, and entropy. On the other hand, 
one should keep in mind that the calculations of these 
quantities at high p may be affected by the fact that the 
measurement of the DOS is limited up to p = 0.572. To 
avoid this effect, we have to set the upper limit as high 
as possible. However, the computation time of the Wang- 
Landau method rapidly increases as the upper limit in- 
creases. 

5. Conclusions 

In this paper, we have presented an efficient Monte 
Carlo method called the LRMC method for lattice glass 
models which are characterized by hard constraint con- 
ditions. Like the TV-fold way method, we make a list of 
the sites into which we can insert a particle, and update 
it whenever the particle configuration is changed. By uti- 
lizing the list, we can avoid a useless trial of an insertion 
which conflicts with the constraint conditions. The use 
of the LRMC method enables us to make the relaxation 
time about 10 3 times shorter than that of the standard 
MC method when p is large. We also have shown that 
the efficiency of the replica exchange method is improved 
much by using the LRMC method as a local update 
method. We can expect such improvement not only in 
the replica exchange method but also in other extended 
ensemble methods. 

In the present work, the LRMC method was applied 
only to the BM model on a regular random graph. How- 
ever, this method is applicable no matter whether the 
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model is defined on a sparse random graph or a usual 
lattice in finite dimensions. An applicable class of mod- 
els includes constraint-satisfaction problems with binary 
variables such as /^-satisfiability problems and vertex 
cover problems. 28 ) We also emphasize that we can use 
the method not only geometrically constrained models 
like the BM model but also kinetically constrained mod- 
els like the Kob- Anderson model. 7 ' 

Another important merit of the LRMC method is that, 
as demonstrated by the scaling plot in Fig. 6(h), the dy- 
namics of the model is not changed essentially by the use 
of the method. We can therefore use the LRMC method 
to investigate the dynamics of the model. It is also worth 
pointing out that this method is applicable even if we 
employ a specific dynamical rule. For example, let us 
consider dynamical rule that a particle can hop from 
a site to one of its nearest-neighbouring sites provided 
that the constraint conditions remain satisfied after the 
hopping. This dynamics is easily feasible in the manner 
of the LRMC method by making and updating the list 
of the nearest-neighbouring pairs where such hopping is 
possible. We hope that this method will stimulate fur- 
ther studies of lattice glass models on both static and 
dynamic properties. 

Acknowledgment 

This work is supported by a Grant-in- Aid for Scientific 
Research (No. 22340109) from the Ministry of Education, 
Culture, Sports, Science and Technology in Japan. 



10 
ii 

12 
13 

14 

15 
16 
17 
18 
1!) 
20 

21 

22 

23 



P. W. Anderson: Science 267 (1995) 1609. 

A. Cavagna: Physics Reports 476 (2009) 51. 
S. P. Das: Rev. Mod. Phys. 76 (2004) 785. 
J. C. Dyre: Rev. Mod. Phys. 78 (2006) 953. 

G. Parisi and F. Zamponi: Rev. Mod. Phys. 82 (2010) 789. 
L. Bcrthier and G. Biroli: Rev. Mod. Phys. 83 (2011) 587. 
W. Kob and H. C. Andersen: Phys. Rev. E 48 (1993) 4364. 
G. Biroli and M. Mezard: Phys. Rev. Lett. 88 (2002) 025501. 
M. P. Ciamarra, M. Tarzia, A. dc Candia, and A. Coniglio: 
Phys. Rev. E 67 (2003) 057105. 

M. Weigt and A. K. Hartmann: Europhys. Lett. 62 (2003) 533. 
C. Toninelli, G. Biroli, and D. S. Fisher: Phys. Rev. Lett. 96 
(2006) 035702. 

M. Mezard and G. Parisi: Eur. Phys. J. B 20 (2001) 217. 
O. Rivoire, G. Biroli, O. C. Martin and M. Mezard: Eur. Phys. 
J. B 37 (2004) 55. 

F. Krzakala, M. Tarzia and L. Zdeborova: Phys. Rev. Lett. 
101 (2008) 165702. 

Y. Iba: 2001 Int. J. Mod. Phys. C 12 (2001) 623. 

B. Berg and T. Neuhaus: Phys. Lett. B 267 (1991) 249. 
B. Berg and T. Neuhaus: Phys. Rev. Lett. 68 (1992) 9. 

F. Wang and D. P. Landau: Phys. Rev. Lett. 86 (2001) 2050. 
F. Wang and D. P. Landau: Phys. Rev. E 64 (2001) 056101. 
K. Hukushima and K. Nemoto: J. Phys. Soc. Jpn. 65 (1996) 
1604. 

A. B. Bortz, M. H. Kalos, and J. L. Lebowitz: J. Comput. Phys. 
17 (1975) 10. 

N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and 
E. Teller: J. Chem. Phys. 21 (1953) 1087. 
We have assumed that the Hamiltonian is a monotone de- 
creasing function of the number of particles. We also have 
assumed that the particle configuration is updated by insert- 
ing or deleting a particle one by one and the particle-hole 
exchange update is not used. On these assumptions, a state 
with K{oi\ = is considered to be in local minimum because 
any movement from the present state, i.e., a deletion of a 



particle from a site, results in an increase in the Hamiltonian. 

24) G. Bhanot, D. Duke, and R. Salvador: Phys. Rev. B 33 (1986) 
7841. 

25) G. Bhanot, D. Duke, and R. Salvador: J. Stat. Phys. 44 (1986) 
985. 

26) K. Hukushima and S. Sasa: J. Phys.: Conf. Scr. 233 (2010) 
012004. 

27) E. Bittner, A. Nufibaumer, and W. Janke: Phys. Rev. Lett. 
101 (2008) 130603. 

28) M. Mezard and A. Montanari: Information, Physics, and 
Computation (Oxford University Press, Oxford, 2009), Chap. 
10, p. 197. 

Appendix A: The method to detect the sites on 
which the insertion list has to be 
updated 

In the LRMC method, we need to update the insertion 
list whenever a particle is inserted into a site or deleted 
from there. In order to do that, we first need to detect 
the sites on which the insertion list has to be updated. In 
this appendix, we describe the method to detect them. 
We emphasize that it is quite important to perform this 
procedure as efficient as possible because this is one of 
the most fundamental procedure in the LRMC method. 

We introduce two variables to explain the method. 
Firstly, a Boolean Insert(p) is .TRUE, if it is possible 
to insert a particle at the site p, or it is .FALSE, oth- 
erwise. The purpose of this appendix is to detect the 
sites at which Insert(p) changes from .TRUE, to .FALSE, 
or conversely by an insertion or deletion of a particle. 
Secondly, an integer NNN(p) denotes the number of the 
nearest-neighbouring particles at site p. Due to the con- 
straint conditions of the BM model, NNN(p) has to be less 
than or equal to I if the site p is occupied by a particle. 
After we choose an initial particle configuration, we set 
Insert (p) and NNN(p) at the beginning of the simulation. 
Then, we update them whenever a particle is inserted or 
deleted. 

The organization of this appendix is as follows: In §A.l, 
we explain basic strategy to detect the sites at which 
Insert(p) changes by an insertion or deletion of a parti- 
cle. In this subsection, we consider the general case that 
the two parameters fc and I of the BM model are arbitrary 
integers. In §A.2, we consider a special case I = 1 and 
explain an optimized method for this case. In these two 
subsections, we assume that the graph does not have a 
loop which involves an insertion or deletion site. We con- 
sider in §A.3 how the method to detect the changed sites 
should be modified when such loop exists. 

A.l Basic strategy in the general BM model 

In this subsection, we consider the general case that 
the two parameters k and I of the BM model are arbi- 
trary integers, and explain the basic strategy to detect 
the sites at which Insert changes from .TRUE, to .FALSE, 
or conversely by an insertion or deletion of a particle. 
Firstly, in the following subsection, we consider at which 
sites such change in Insert may occur. It is important 
to reduce the possible sites as many as possible to reduce 
the computational time. Secondly, in §A.1.2, we explain 
the general method to judge whether a change in Insert 
really occurs or not. 
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Fig. A-l. (Color online) An example when Insert (r) depends on 
(T p , where r is a next nearest-neighbouring site of p. The values 
of the two parameters of the BM model are k = 3 and 1 = 1. In 
the particle configuration (b), it is impossible to insert a particle 
into the site r because the number of the nearest-neighbouring 
particles of the site q becomes larger than Z = 1 by the insertion. 

A. 1.1 Specification of the possible sites 

We consider to insert or delete a particle at a site p. 
Then, g p changes from to 1 or conversely. There are 
three kinds of sites where Insert may change: the site 
p itself, the nearest- neighbouring sites of p, and the next 
nearest-neighbouring sites of p. Firstly, let us consider 
the case of the site p itself. After wc delete a parti- 
cle from the site p, Insert(p) becomes .TRUE.. On the 
other hand, after we insert a particle into the site p, 
Insert(p) becomes .FALSE, because the site p has al- 
ready been occupied by a particle. These two facts mean 
that Insert (p) always changes by the insertion or dele- 
tion of a particle at the site p. Therefore, we do not need 
to judge whether the change occurs or not. Secondly, 
concerning a nearest-neighbouring site q, it is apparent 
that Insert (q) is .FALSE, if the site q is already occupied 
by a particle. This means that we do not need to judge 
whether Insert(g) changes or not if o~ q is 1. Lastly, we 
consider the case of a next nearest-neighbouring site r. 
As it is shown in Fig. A-l, if Insert(r) depends on <r p , 
the nearest-neighbouring site q which is between p and r 
has to be occupied by a particle. In other words, we do 
not need to judge whether Insert(r) changes or not if 
o q is 0. 

By taking all of these considerations into account, we 
find that the following procedure is enough to list all 
of the sites where Insert changes by the insertion or 
deletion of a particle at the site p: 

(I) Add the site p into the list. 

(II) Perform the following procedure for all of the 
nearest-neighbouring sites q: 

(i) If a q = 0, judge whether Insert(g) changes or 
not. If it changes, add the site q into the list. 

(ii) If G q = 1, check all of the next nearest- 
neighbouring sites r which connect with the site q 
whether Insert(r) changes or not. If it changes, 
add the site r into the list. 

A. 1.2 General method to judge whether Insert changes 
or not 

When Insert(g) may change at a site q by the inser- 
tion or deletion of a particle at the site p, we have to judge 
whether the change really occurs or not. In this sub- 
section, we consider how we judge it. Because, as men- 
tioned above, Insert is updated whenever a particle is 
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Fig. A-2. (Color online) Three examples when one of the three 
necessary and sufficient conditions for Insert (i) = .TRUE, is not 
satisfied, and two examples when all of the three conditions are 
satisfied. The values of the two parameters of the BM model are 
k = 3 and 1 = 1. The first, the second, and the third condition 
are not satisfied in the case (i), (ii), and (iii), respectively. In 
contrast, all of the three conditions are satisfied in the cases (iv) 
and (v). 

inserted or deleted, we know whether Insert (q) is .TRUE, 
or .FALSE, before the insertion or deletion. Therefore, to 
judge whether the change occurs or not, it is sufficient 
to know whether Insert(g) is .TRUE, or .FALSE, after the 
insertion or deletion. As shown in Fig. A-2, the neces- 
sary and sufficient conditions for Insert(i) = .TRUE, arc 
as follows: 

• Gi = 0. 

• NNN(i) < I. 

• NNN(j) < I — 1 for all of the nearest-neighbouring 
sites j of the site i which are occupied by a particle. 

Therefore, Insert (i) is .TRUE, if all of the three condi- 
tions are satisfied. Otherwise, it is .FALSE.. 

A. 2 Optimized method for I = 1 

In principle, we can judge whether Insert changes 
or not by the method described in §A.1.2. However, it 
is time-consuming to examine whether all of the three 
necessary and sufficient conditions mentioned in §A.1.2 
arc satisfied or not for all of the sites at which Insert 
may change. In this subsection, we consider a special 
case 1 = 1, and consider to reduce the procedure for the 
judgement as much as possible. 

Before we describe the details of the method, for the 
convenience of explanation, we show the necessary and 
sufficient conditions for Insert(i) = .TRUE, in the case 
of 1 = 1: 

(a) Gi = 0. 

(b) NNN(i) is either or 1. 

(c) If NNN(i) is 1, NNN(j) = 0, where j is the nearest- 
neighbouring site which is occupied by a particle. 

A. 2.1 Judgement at a next nearest- neighbouring site 

In this subsection, we consider to judge whether 
Insert(r) changes or not at a next nearest-neighbouring 
site r. Figure A-3 shows an example of possible parti- 
cle configurations when we judge the site r. The particle 
configuration changes from (a) to (b) or conversely by 
inserting or deleting a particle at the site p, respectively. 
We see that these particle configurations satisfy the fol- 
lowing three conditions: 

(1) G q = 1, where q is the nearest- neighbouring site be- 
tween p and r. 
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Fig. A-3. (Color online) An example of particle configurations 
when we judge next nearest-neighbouring sites. The values of 
the two parameters of the BM model are k = 3 and 1 = 1. The 
particle configuration changes from (a) to (b) or conversely by in- 
serting or deleting a particle at the site p, respectively. Insert (r) 
changes by the insertion or deletion of a particle at the site p, 
whereas such change does not occur at the site s. 

(2) Both oy and o~ s is 0, where r and s are the next 
nearest-neighbouring sites which connect with q. 

(3) Both Insert(r) and Insert(s) are .FALSE, when 
(Tp = 1. 

The first condition comes from the fact that we need to 
judge the site r only if a q is 1 (see the procedure (ii) 
in the last paragraph of §A.1.1). The second condition 
comes from the fact that the constraint condition has 
to be satisfied at the site q even if a p = 1. The third 
condition comes from the fact that NNN(g) has already 
been 1(= I) when a. p = 1. From the condition (3), we 
find that Insert(r) changes by inserting or deleting a 
particle at the site p if and only if Insert(r) is .TRUE, 
when Op = 0. 

Therefore, we next consider the necessary and suffi- 
cient conditions so that Insert(r) is .TRUE, under the 
prior condition a p = 0. As we mentioned above, there are 
the three necessary and sufficient conditions (a)-(c) men- 
tioned at the beginning of of §A.2 so that Insert(r) is 
.TRUE.. However, in practice, we do not need to check all 
of the three conditions. Firstly, the condition (a) is satis- 
fied automatically due to the condition (2). Secondly, if 
NNN(r) is 1 and the condition (b) is satisfied (note that 
it is impossible that NNN(r) = because a q = 1), the 
condition (c), i.e., NNN(g) = 0, is also satisfied automati- 
cally because a p = by the prior condition and both ay 
and a s are by the condition (2). From these two facts, 
we find that the condition (b) is the only necessary and 
sufficient condition we have to check. 

In conclusion, Insert (r) changes by the insertion or 
deletion of a particle at the site p if and only if NNN(r) is 
1. 

A. 2. 2 Judgement at a nearest-neighbouring site after 
the insertion of a particle into the site p 
In this subsection, we consider to judge whether 
Insert(q) changes or not by the insertion of a particle at 
the site p, where q is a nearest-neighbouring site. Because 
we insert a particle into the site p, it is impossible that 
Insert(g) changes from .FALSE, to .TRUE, by the inser- 
tion. Therefore, the necessary and sufficient conditions 
that Insert (q) changes are that Insert (q) before the in- 
sertion is .TRUE, and that Insert (q) after the insertion 
is .FALSE.. Because Insert is updated at each step, we 



know Insert(q) before the insertion. Therefore, we can 
easily check whether the first condition is satisfied or 
not. We next consider the second condition. To this end, 
it is convenient to consider the necessary and sufficient 
conditions for the complementary event that Insert(g) 
after the insertion is .TRUE.. In general, the necessary 
and sufficient conditions are the three conditions (a)-(c) 
mentioned at the beginning of §A.2. However, because 
we need to judge the site q only if a q is (see the pro- 
cedure (i) in the last paragraph of §A.1.1), the condition 
(a) is always satisfied. We also notice that NNN(g) can not 
be after the insertion at the site p. Therefore, the nec- 
essary and sufficient conditions that Insert(g) after the 
insertion is .TRUE, arc that NNN(q) = 1 and NNN(p) = 0. 
This means that Insert (q) after the insertion is .FALSE, 
if and only if either NNN(g) ^ 1 or NNN(p) ^ 0, or both. 

In conclusion, Insert(q) changes from .TRUE, to 
.FALSE, if and only if the following two conditions are 
satisfied: 

• Insert(g) is .TRUE, before a particle is inserted at 
the site p. 

• Either NNN(q) ^ 1 or NNN(p) ^ 0, or both, after a 
particle is inserted into the site p. 

It should be noted that NNN(g) and NNN(p) are the num- 
bers of the nearest-neighbouring particles after a particle 
is inserted into the site p. 

A. 2.3 Judgement at a nearest-neighbouring site after 
the deletion of a particle from the site p 

In this subsection, we consider to judge whether 
Insert(g) changes or not by the deletion of a particle at 
the site p, where q is a nearest-neighbouring site. From a 
consideration similar to that in the previous subsection, 
we find that the necessary and sufficient conditions that 
Insert(g) changes are that Insert^) before the dele- 
tion is .FALSE, and that Insert^) after the deletion is 
.TRUE.. Firstly, as mentioned in the previous subsection, 
we can easily check whether the first condition is satisfied 
or not. Secondly, Insert(q) after the deletion is .TRUE, if 
and only if all of the three conditions (a)-(c) is satisfied. 
However, since the site q is judged only if a q is (see the 
procedure (i) in the last paragraph of §A.1.1), we can 
remove the condition (a) from them. 

In conclusion, Insert (q) changes from .FALSE, to 
.TRUE, if and only if the following three conditions are 
satisfied: 

• Insert(g) is .FALSE, before a particle is deleted from 
the site p. 

• NNN(g) is or 1. 

• If NNN(g) is 1, NNN(j) is 0, where j is the nearest- 
neighbouring site of q which is occupied by a parti- 
cle. 

We again emphasize that NNN(q) and NNN(j) are the num- 
bers of the nearest-neighbouring particles after a particle 
is deleted from the site p. 

A. 3 The case when a loop exists around the insertion 
or deletion site 
In §A.l and §A.2, we have implicitly assumed that 
there is no loop which involves the insertion or dele- 
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Fig. A-4. (Color online) Three examples of a graph when a loop 
which involves the insertion or deletion site p exists. The length 
of the loop is 3 in the case (a), 4 in the case (b), and 5 in the 
case (c). 



tion site p. In this subsection, we explain how we should 
change the method described in §A.l and §A.2 when a 
loop which involves the site p exists. Figure A-4 shows 
examples of a graph when there is a loop. The length 
of the loop is 3 in the case (a), 4 in the case (b), and 5 
in the case (c). We can check that the method to judge 
whether Insert changes or not is still valid even if such 
loop exists. Therefore, we do not need to modify both 
the general method described in §A.1.2 and the opti- 
mized method for I = 1 described in §A.2. The only part 
we should change is the specification of the possible sites 
described in §A.1.1. 

When a loop whose length is either 3 or 4 exists, i.e., in 
the cases (a) and (b) in Fig. A-4, we should be careful not 
to check a site twice as a possible site. When the length 
of the loop is equal to or larger than 5 like the case (c), 
we do not need to care about it. We hereafter consider 
the two cases (a) and (b). Firstly, when we check the site 
q in the case (a), we should exclude the site r from the 
next nearest-neighbouring sites so as not to check this 
site in the procedure (ii) in the last paragraph of §A.1.1. 
Otherwise, the site r is checked twice because this site 
is checked in the procedure (i) as a nearest-neighbouring 
site. To avoid this double check, in the procedure (ii), we 
should check a site r when it satisfies the following two 
conditions: 

• r is a nearest-neighbouring site of q. 

• r is neither p nor one of the nearest-neighbouring 
sites of p. 

Because these two conditions only depend on the shape 
of the graph, it is enough to calculate the set of sites 
which satisfy the two conditions once at the beginning 
of the simulation. This calculation should be done for 
each site. We next consider the case (b) in Fig. A-4. We 
suppose that, as shown in the figure, both q and r are 
occupied by a particle. Then, if we naively perform the 
procedure described in §A.1.1, the site t is checked twice 
as a possible site. When 1 = 1, this kind of double check 
does not occur even if we naively perform the procedure 
because it is impossible that both q and r are occupied by 
a particle. Note that NNN(p) < I = 1 because Insert(p) is 
.TRUE, when a p = 0. However, if I > 2, we should modify 
the procedure so as to check the site t once when both q 
and r are occupied by a particle. 



Appendix B: The method to update the inser- 
tion list 

In this appendix, we explain the method to update 
the insertion list. Before we start the explanation, we 
explain the situation when we update the insertion list 
and introduce several technical terms. We assume that 
we know all of the sites into which we can insert a par- 
ticle. This information is stored in an array List. The 
value of List(m) denotes the m-th site into which we 
can insert a particle. The number of the insertable sites 
is stored in an integer Nlist. We assume that List(m) is 
for m > Nlist. We also assume that we know where an 
insertable site i is recorded in the array List. This infor- 
mation is stored in an array Reverse(t) and the value of 
Reverse(i) denotes the position in the array List. That 
is to say, if the number of Reverse(i) is, say, 10, the site i 
is insertable and List(10) = i. It is worth noticing that, 
if we set Reverse(i) to a negative integer such as — 1 
when the site i is not insertable, we do not need to pre- 
pare the array Insert(i) introduced in the appendix A 
because we can store the information whether the site i 
is insertable or not in Reverse(i). In this setting, posi- 
tive Reverse(i) corresponds to Insert(z) = .TRUE, and 
negative Reverse(i) corresponds to Insert(i) = .FALSE.. 
After we choose an initial particle configuration, we set 
the two arrays List(m) and Reverse(i) and the inte- 
ger Nlist at the beginning of the simulation. Then, we 
update them whenever a particle is inserted or deleted. 

We now start to explain the method to update the 
insertion list. We assume that we have detected all of 
the sites on which the insertion list has to be updated by 
using the method described in the appendix A. There 
are two operations to update the insertion list: 

(A) Remove a site i from the list. 

(B) Add a site i into the list. 

We first consider the operation (A). Now the point is 
that, when we insert a particle, we choose the insertion 
site at random from the array List. This means that the 
order in the list is not important. Therefore, when we 
remove the m-th element from the list (m = Reverse(i)), 
we can replace the m-th element with the last element 
and set the last element to zero after that. By taking this 
and the fact that Reverse(i) = m into account, we can 
perform the operation (A) in the following way: 

Reverse (List (Nlist) )=Reverse(i) 

List (Reverse (i) ) =List (Nlist) 

List(Nlist)=0 

Nlist=Nlist-l 

Reverse(i)=-1 

The operation (B) is simpler. We just add the site i at 
the end of the list. This is performed in the following 

way: 

Nlist=Nlist+l 
List(Nlist)=i 
Reverse (i)=Nlist 



